***Remember the order of variables is X, Z, W. ***
***On the graph, W represents the different lines.  I set them to W0=-1SD PAGG W1=50% W2=+1SD W3=95%***

* two-sided, p<.05=1.96, p<.1=1.645
* y = x z w xz xw zw xzw*
* x(motive), z(aggress), w(Viol), y(dependent variable)
*Estimating marginal effect of Violent Treatment, conditional on TA and group attitude
reg vote viol aggress commitment violxaggress violxcommitment aggressxcommitment violxaggxcommit

#delimit ;
*Estimating marginal effect of group attitude, conditional on TA and treatment*;

egen zmin = min(aggress);
egen zmax = max(aggress);

gen z0 = (((_n-1)/(20-1))*(zmax-zmin))+zmin in 1/20;
generate Z=z0;

*     ****************************************************************  *;
*       Generate the values of W for which you want to calculate the    *;
*       marginal effect (and standard errors) of X on Y.                *;
*     ****************************************************************  *;

scalar W0=0;
scalar W1=1;

*     ****************************************************************  *;
*       Grab elements of the coefficient and variance-covariance matrix *;
*       that are required to calculate the marginal effect and standard *;
*       errors.                                                         *;
*     ****************************************************************  *;

matrix b=e(b);
matrix V=e(V);

scalar b1=b[1,1];
scalar b2=b[1,2];
scalar b3=b[1,3];
scalar b4=b[1,4];
scalar b5=b[1,5];
scalar b6=b[1,6];
scalar b7=b[1,7];

scalar varb1=V[1,1];
scalar varb2=V[2,2];
scalar varb3=V[3,3];
scalar varb4=V[4,4];
scalar varb5=V[5,5];
scalar varb6=V[6,6];
scalar varb7=V[7,7];

scalar covb1b4=V[1,4];
scalar covb1b5=V[1,5];
scalar covb1b7=V[1,7];
scalar covb4b5=V[4,5];
scalar covb4b7=V[4,7];
scalar covb5b7=V[5,7];

scalar list b1 b2 b3 b4 b5 b6 b7 varb1 varb2 varb3 varb4 varb5 varb6 varb7 
           covb1b4 covb1b5 covb1b7 covb4b5 covb4b7 covb5b7;
           
*     ****************************************************************  *;
*       We want to calculate the marginal effect of X on Y for all      *;
*       Z values of the modifying variable Z. We also want to          *;
*       calculate this marginal effect as Z changes for specific values *;
*       of the second modifying variable W.  In the code below, we      *;
*       calculate the marginal effect of X on Y for all values of Z     *;
*       when W=0, when W=1, when W=2, and when W=3.                     *;
*     ****************************************************************  *;

gen conb0=b1+b4*Z+b5*W0+b7*(Z*W0);
gen conb1=b1+b4*Z+b5*W1+b7*(Z*W1);

*     ****************************************************************  *;
*       Calculate the standard errors for the marginal effect of X on Y *;
*       for all Z values of the modifying variable Z. Do this for the  *;
*       case when W=0, when W=1, when W=2, and when W=3.                *;
*     ****************************************************************  *;

gen conse0=sqrt(varb1
               + varb4*(Z^2) + varb5*(W0^2) + varb7*(Z^2)*(W0^2)
               + 2*Z*covb1b4 + 2*W0*covb1b5 + 2*Z*W0*covb1b7 + 2*Z*W0*covb4b5
               + 2*W0*(Z^2)*covb4b7 + 2*(W0^2)*Z*covb5b7);

gen conse1=sqrt(varb1
               + varb4*(Z^2) + varb5*(W1^2) + varb7*(Z^2)*(W1^2)
               + 2*Z*covb1b4 + 2*W1*covb1b5 + 2*Z*W1*covb1b7 + 2*Z*W1*covb4b5
               + 2*W1*(Z^2)*covb4b7 + 2*(W1^2)*Z*covb5b7);
                
			   
*     ****************************************************************  *;
*                           Create t statistics                         *;
*     ****************************************************************  *;

gen t0=conb0/conse0;
gen t1=conb1/conse1;

*     ****************************************************************  *;
*       Generate a `shadow' variable that is equal to the marginal      *;
*       marginal effects that have already been calculated.  These will *;
*       be helpful in placing the stars in the figure as you will see.  *;
*     ****************************************************************  *;

gen consb0=conb0;
gen consb1=conb1;

gen consb090=conb0;
gen consb190=conb1;

*     ****************************************************************  *;
*       Make the `shadow' variable be missing if the t score is not     *;
*       larger than the critical level of significance that you want.   *;
*     ****************************************************************  *;

replace consb0 = . if abs(t0)<1.645;
replace consb1 = . if abs(t1)<1.645;

replace consb090 = . if abs(t0)<1.28|abs(t0)>=1.645;
replace consb190 = . if abs(t1)<1.28|abs(t1)>=1.645;

*     ****************************************************************  *;
*       Generate a string variable called str1 that is designated with  *;
*       a star.                                                         *;
*     ****************************************************************  *;

generate str1 txt="*";
generate str2 txt2="+";

*     ****************************************************************  *;
*       Graph the marginal effect of X on Y across the desired range of *;
*       the modifying variable Z. Do this for when W=0 and for when W=1.*;
*                                                                       *;
*     ****************************************************************  *;


graph twoway   line conb0  Z,  lpattern(dash) lwidth(medium)
       ||  scatter consb0 Z,  mlabel(txt) msymbol(i) mlabsize(4) mlabgap(-1.0) mlabposition(11)
       ||     line conb1  Z,  lpattern(solid) lwidth(medium)
       ||  scatter consb1 Z,  mlabel(txt) msymbol(i) mlabsize(4) mlabgap(-1.0) mlabposition(11)  
	   ||  scatter consb090 Z,  mlabel(txt2) msymbol(i) mlabsize(4) mlabgap(-1.0) mlabposition(11)
       ||  scatter consb190 Z,  mlabel(txt2) msymbol(i) mlabsize(4) mlabgap(-1.0) mlabposition(11)  
       ||  ,   
           xlabel(, labsize(5)) 
           ylabel(, labsize(5))
           yscale(noline)
           xscale(noline)
 xsize(3)
           legend(order(1 "Low Motivation" 2 "High Motivation"))
           xtitle(Trait Aggression, size(5))
           ytitle(Marginal Treatment FX on Predicted Turnout Values, size(5))
           xsca(titlegap(2)) 
           ysca(titlegap(4))
           scheme(s2mono) graphregion(fcolor(white));
       
	   
drop zmax zmin z0 Z conb0 conb1 consb090 consb190 conse0 conse1  t0 t1 consb0 consb1 txt txt2;
